Spatiotemporal evolution of polaronic states in finite quantum systems 
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We study the quantum dynamics of small polaron formation and polaron transport through finite 
quantum structures in the framework of the one-dimensional Holstein model with site-dependent 
potentials and interactions. Combining Lanczos diagonalization with Chebyshev moment expansion 
of the time evolution operator, we determine how different initial states, representing stationary 
ground states or injected wave packets, after an electron-phonon interaction quench, develop in 
real space and time. Thereby, the full quantum nature and dynamics of electrons and phonons 
is preserved. We find that the decay out of the initial state sensitively depends on the energy 
and momentum of the incoming particle, the electron-phonon coupling strength, and the phonon 
frequency, whereupon bound polaron- phonon excited states may emerge in the strong-coupling 
regime. The tunneling of a Holstein polaron through a quantum wall or dot is generally accompanied 
by strong phonon number fluctuations due to phonon emission and reabsorption processes. 



I. INTRODUCTION 

Electrons injected into low-dimensional quantum 
structures with strong electron-phonon (EP) interac- 
tion can cause local lattice deformations and thereby 
relax to "self-trapped" polaron states. Polaron self- 
trapping does not imply a breaking of translational in- 
variance. When the polaron forms, the electron is be- 
ing dressed by a phonon cloud, and the collective state 
translates through the lattice. This indicates that vi- 
brational modes for polaron transport through nanoscale 
quantum devices are of vital importance. The micro- 
scopic structure of polarons and the contexts in which 
they appear are rather diverse. Phonon and polaron ef- 
fects have been investigated for , e.g., molecular transis- 
tors^ quantum dots^ tunneling diodes and Aharonov- 
Bohm rings, 3 metal/organic/metal structures^ and Car- 
bon nanotubes, 5,6 although primarily with respect to 
steady state properties. Recently, time-resolved spec- 
troscopy has made it possible to address also the dy- 
namical aspects of self-trapping, e.g., by directly time re- 
solving the vibrational motions associated with the local- 
ized carrier. Taking advantage of the ultra-short pulse- 
widths of recent lasers, the femto-second dynamics of po- 
laron formation and exciton-phonon dressing has been 
observed in pump-probe experiments^ 

From a theoretical point of view, describing the time 
dependence of small polaron formation requires a physics 
that is related to particle and phonon dynamics on the 
scale of the unit celli* The simplest model that captures 
such a situation is the lattice-polaron Holstein model. 9 
This model assumes that the orbital states are identical 
on each site and the particle can move from site to site 
exactly as in a tight-binding model. The phonons are 
coupled to the particle at whichever site it is on. The 
dynamics of the phonons is treated purely locally with 
Einstein oscillators representing the intra-site (molecu- 
lar) vibrations. After six decades of intense research the 
equilibrium properties of the Holstein model are well- 
understood, at least in the single-particle sector (for re- 



cent reviews, see Refs. [TolfTll ). In contrast, there is a 
rather incomplete understanding of time-dependent and 
out-of-equilibrium phenomena. Some issues seem to be 
settled, e.g., the charge-transfer and correlated charge- 
deformation dynamics (but only for a two-site Holstein 
model) jia the increase of the polaron formation time in 
two and three dimensions due to an adiabatic potential 
barrier between extended electron and self-trapped po- 
laron states; 13 ! 14 or the conditional hopping rate of an 
injected electron and the vibrational relaxation time^ 
Quite recently the problem of determining the polaron 
formation time has been tackled J£ Many questions re- 
main at least partly unsolved however. For instance, in 
what way is quantum-dot polaron formation different in 
the adiabatic and anti-adiabatic regimes? And how does 
a bare particle evolves into a polaron after an "interac- 
tion quench"? Or, how does a polaronic quasiparticle 
tunnel through a potential barrier? 

In this paper, we address some of these questions. 
To that end, we calculate — by means of numerically 
exact Lanczos diagonalization and Chebyshev expan- 
sion techniques — the real space and time evolution of 
polaronic states in the one-dimensional Holstein model 
with spatially and temporally varying on-site potentials 
and/or EP interaction strengths. The proposed approach 
is applicable to the dynamics of quasiparticle formation 
in several branches of physics. 



II. MODEL AND METHOD 



Modified Holstein Hamiltonian 



Focusing on polaron formation in one-dimensional fi- 
nite quantum structures with short-range non-polar EP 
interaction we consider the generalized Holstein molecu- 



2 



lar crystal model^^ 

i i 



(1) 



where cj (cj and b\ (6J are creation (annihilation) op- 
erators for electrons and dispersionless optical phonons 
on site i, respectively, and n i = c]c^ is the corresponding 
particle number operator. In (TTJ), the site-dependent po- 
tentials Aj can describe a tunnel barrier, a voltage bias, 
or disorder effects, to denotes the nearest-neighbor elec- 
tron transfer integral, and gi gives the local interaction 
of an electron on Wannier site i to an internal vibrational 
mode with frequency wo- 

The ratio oJo/tg determines which of the two subsys- 
tems, electrons or phonons, is the fast or the slow one. 
In the adiabatic limit (wo/io) "C 1, the motion of the 
particle is affected by quasi-static lattice deformations, 
whereas in the opposite, anti-adiabatic limit (woAo) 3> 1 
the lattice deformation is presumed to adjust instanta- 
neously to the position of the carrier. 

The dimensionless EP coupling constant g 2 normally 
appears in (small polaron) strong-coupling perturbation 
theory, where it describes the polaronic mass enhance- 
ment m*/m = e 9 (for homogeneous systems, gi = g). 
There is another natural measure of the strength of the 
electron-phonon interaction, the familiar polaronic level 
shift E p . At strong EP coupling, E p gives the leading- 
order energy shift of the band dispersion. 18 In general, 
there is no simple relation between g 2 and E p . If the EP 
coupling is local and the phonon mode is dispersionless, 
however, then g 2 = E p /loq, and E p is usually identified 
with the polaron binding energy*^ 

The crossover from essentially free electronic carriers 
to heavy polaronic quasiparticles is known to occur for a 
translational invariant system, provided that two condi- 
tions, g 2 > 1 and E p /zto > 1 (z = 2 ( in one dimension), 
are fulfilled iSS So while the first requirement is more re- 
strictive in the anti-adiabatic case, the formation of a 
small polaron state will be determined by the second cri- 
terion in the adiabatic regime. This likewise holds for 
the generalized Holstein model JTJ where g 2 = E Pyi /uJo 
[with a view to the different cases studied in Sec. IIII 
we split E Pt i — e p + e Pt i up into a constant (e p ) and a 
site-dependent part 

When investigating the physically most interesting 
crossover regime of the Holstein model where polarons 
form, i.e. the self-trapping transition of the charge car- 
riers takes place, standard analytical approaches fail to 
a large extent. This is because, precisely in this situ- 
ation, the characteristic electronic and phononic energy 
scales are not well separated. So far quasi approximation- 
free numerical methods like quantum Monte Carlo sim- 
ulations^ exact diagonalizations, 22 or density-matrix 
renormalisation group techniques 2 ^ yield the most reli- 
able results for the ground-state and spectral properties 



of Holstein polarons. 



B. Chebyshev expansion technique 

To study the real space and time formation of a 
polaronic quasiparticle from a bare electron the time- 
dependent many-body Schrodinger equation has to be 
solved. For systems with moderate Hilbert space di- 
mensions a full diagonalization of the Hamiltonian al- 
lows for an exact calculation of the quantum state at 
arbitrary times. Because of the phonon degrees of free- 
dom the Hilbert space of the Holstein model is infinite, 
even for a finite lattice and in the single-particle sec- 
tor. Truncating the Hilbert space of the phonons or con- 
structing a variational Hilbert space including multiplc- 
phonon excitations ; 10 ' 24 a direct numerical integration of 
the Schrodinger equation can be performed, yielding the 
polaron many-body wave function at early times^ Al- 
ternatively one can exploit a Chebyshev moment based 
expansion of the time evolution operator*^ Since this 
technique also applies to very general situations and has 
been proven to be superior to direct integration and other 
iterative Schrddinger-equation solution schemes as to its 
efficiency (i.e. computational costs) and accuracy^ the 
remainder of this section briefly outlines this less well- 
known approach. 

The time evolution of a quantum state \ip) is described 
by the Schrodinger equation 



ih^\m)=mm) 



(2) 



If the Hamilton operator H does not explicitly depend 
on time t we can formally integrate this equation and 
express the dynamics of an initial state |^(0)) in terms 
of the time evolution operator U (t, 0) as 



\il>(t)) = U(t,0)MQ)), 



where 



J7(t,0) 



iHt/n 



(3) 



(4) 



The time evolution operator U(t + At, t) = U (At) for 
a given (usually small) time step At can be expanded in 
a finite series of Nc first-kind Chebyshev polynomials of 
order n, 



T n (x) = cos[narccos(x)] 



(5) 



We obtainS§r^ 
U{At) = e - ibA */ ri 



CO 



(aAt/h) + c n {aAt/K)T n (H) 



n=l 



(6) 

Prior to the expansion, the Hamiltonian has to be 
shifted and rescaled such that the spectrum of H — 
(H — b)/a is within the definition interval of the Cheby- 
shev polynomials, [—1,1]^ The parameters a and b 
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are calculated from the extremal eigenvalues of H as 
b = |(£ m ax + E min ) and a = ±(E max - E min + e). Here 
we introduced e = a(E max — E m - ln ) to ensure the rescaled 
spectrum \E\ < l/(l + a) lies well inside [—1, 1]. In prac- 
tice, we use a = 0.01. The Chebyshev expansion also 
applies to systems with Holstein-type unbounded spec- 
tra.— Here we can truncate the infinite Hilbert space to 
a finite dimension by restricting the model on a discrete 
space grid or using an energy cutoff. In this way we en- 
sure the finiteness of the extreme eigenvalues. 
In ([5]), the expansion coefficients c n are given by 

" l ' = dx = (-i) n J n (aAt/h) ; 
7rvl — x 2 

(7) 

J n denotes the ri-th order Bcsscl function of the first kind. 

In order to calculate the evolution of a state \i/j(t)) from 
one time grid point to the adjacent one, 

\iP(t + At))=U(At)\iP(t)), (8) 
we have to accumulate the c„-weighted vectors 

\v n )=T n (H)\m)- (9) 

Since the coefficients c n (aAt/h) depend on the time step 
but not on time explicitly, we need to calculate them 
only once. Instead of evaluating Eq. ([5]) with x — H, the 
vectors \v n ) can be computed iteratively exploiting the 
recurrence relation of the Chebyshev polynomials, 

\v n+1 ) = 2H\v n ) - K_x) , (10) 

with \vi) = H\vq) and \vq) = \ip(t)). Evolving the wave 
function from one time step to the next requires Nq ma- 
trix vector multiplications (MVMs) of a given complex 
vector with the sparse Hamilton matrix of dimension D . 
Of course, to proceed from t = to t, the procedure has 
to be performed t/ At times. 

Note that such a Chebyshev expansion may also be ap- 
plied to systems with time-dependent Hamiltonians, but 
there the time variation of H(t) determines the maxi- 
mum At by which the system may be propagated in a 
single time step. For time-independent H, in principle, 
arbitrary large time steps are possible at the expense of 
increasing Nq ■ We may choose Nq such that for n > Nq 
the modulus of all expansion coefficients 

\c n (aAt/h)\~ J n (aAt/h) (11) 

is smaller than a desired accuracy cutoff. This is facil- 
itated by the fast asymptotic decay of the Bessel func- 
tions, 

1 / GaAt \ n 

JJaAb/h) == [ ^— ) for n ->• oo . (12) 

V27m V 2nn / 

Hence for 2hn 3> eaAt the expansion coefficients c„ de- 
cay superexponential and the series can be truncated 



with negligible error^ In the numerics of Sec. Ill, we 
work with Nq > 10, such that the last moment retained 
\Jn c \ < 10 -9 , i.e. the Chebyshev expansion can be con- 
sidered as quasi-exact, and permits a considerably larger 
time step than e.g. the Crank-Nicholson schemei 26 ' 30 Of 
course, the ground-state energy Eo(t) is unaltered during 
the simulation time. 

Besides the high accuracy of the method, the linear 
scaling of computation time with both time step and 
Hilbert space dimension are promising in view of poten- 
tial applications to more complex systems. Here almost 
all computation time is spent in sparse MVMs, which can 
be efficiently parallelized, allowing for a good speedup on 
parallel computers. We use a memory saving implemen- 
tation of the MVM where the non-zero matrix elements 
are not stored but recomputed in each sparse MVM step, 
limiting the overall memory consumption of our imple- 
mentation to five vectors of size D. In this context we 
can access a massively parallel sparse MVM code which 
has proven to be sufficient to compute the ground state 
of the model (JlJ up to D = 3.5 x 10 11 very efficiently on 
more than 5000 processor cores^i For the single polaron 
dynamics presented here, the matrix dimension is about 
D = 6 x 10 8 and we run the Chebyshev approach on 18 
processors of an SGI Altix4700 compute server, access- 
ing a total of approximately 60 GB of main memory and 
consuming less than 1500 CPU-hrs to compute e.g. the 
results presented in Sec. III. C. 



III. NUMERICAL RESULTS AND DISCUSSION 

In this section we combine exact diagonalization and 
Chebyshev expansion methods; 24 ' 25 ! 29 working in the 
tensorial product Hilbert space of electrons and phonons. 
We set h = 1 and give all energies in units of to. The 
time t will be measured with respect to the character- 
istic electronic and phononic time scales r e = t^ 1 and 
r p h = (wo/27r) _1 , respectively. We consider the case of a 
single electron only. 

A. Interaction quench 

To understand the basic features of the polaron for- 
mation process in the time domain, we first study a 
single oscillatory site to which the Holstein molecular 
crystal model applies, sandwiched between two "wires" 
where electrons are not coupled to phonons (e p = 0). 
The system size is N — 17 with open boundary condi- 
tions (OBCs) at sites i = 1, 17. We consider the case 
Aj = 0. The deformable site is located midway, i = 9. 
Before time t = the system is assumed to be in the 
non-interacting (free-electron) ground state; its energy is 
E a = -1.9696. Then, at t = the EP interaction at 
site 9 is abruptly switched to a positive value e Pj g (an 
interaction "quench" ) , what means that the electron and 
phonon subsystems are locally linked hereafter. Since the 
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whole system is isolated from the environment the total 
energy is conserved during the quench. 

The time evolution of various quantities after such an 
interaction quench is shown in Figs. [1] to [9] for charac- 
teristic situations, ranging from weak to strong EP cou- 
pling and adiabatic to anti-adiabatic cases. As can be 
seen from Figs. [T] to [3] the quantum dynamics after the 
quench depends on the EP coupling strength and phonon 
frequency in a very sensitive way. 



Adiabatic regime 




FIG. 1: (Color online) Time dependence of the particle den- 
sity (upper panel) and phonon number (lower panel) at the 
deformable "molecular crystal" site 9, starting out from the 
free-electron ground state of a 17 site chain with OBCs. The 
phonon frequency uio = 0.3. The EP coupling e p ,9 is switched 
on at t — 0, different curves belong to e p ,g values increased 
by 0.1. In the numerical calculations we take into account up 
to M = 100 phonons, TVc = 10 Chebyshev moments, and use 
a time step At/r e = 0.01. 



Figure Q] illustrates the time evolution of the particle 
density at the oscillatory site after the interaction quench 




FIG. 2: (Color online) Time dependence of the particle den- 
sity (upper panel) and phonon number (lower panel) at molec- 
ular crystal site 9 in the intermediate EP coupling adiabatic 
regime. The initial state is the same as in Fig. [T] e Pt g now is 
increased in steps of 0.04. 



for weak-to-intermediate EP couplings and phonon fre- 
quencies wo smaller than the electronic transfer integral 
to. We note that the electron is not uniformly spread 
over the lattice even at t < where e Pi g = because of 
of the OBCs: (ng)(0) = 1/9 is roughly twice the mean 
electron density 1/17. 

After the local EP interaction is turned on the electron 
can couple to the molecular vibrations at site 9. The ba- 
sic interaction process is the absorption and emission of 
a phonon by the electron with a simultaneously change 
of the electron state. At the same time the lattice is 
distorted locally. Such a lattice distortion may trap the 
charge carrier if the EP coupling is strong. As a result 
the local particle density is enhanced. Since the trap- 
ping potential itself depends on the carrier's state, this 
highly non-linear feedback phenomenon is called "self- 



trapping 



» 32.33 



Figure [T] clearly shows an initial strong 



increase of the local particle density in time. Since the 
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characteristic nearest-neighbor hopping time of a bare 
electron is r e , all "electrons" initially moving toward the 
central site will reach this site within t/r e < 8 (deal- 
ing with a single particle we actually thereby think of 
electronic contributions). This explains the small hump 
on the left shoulder of the first (ng)-increase at about 
t/r e ~ 8. Electrons that move away from the central 
site will reach it after reflection at the boundaries within 
the time interval 8 < t/r e < 17. So if the particle 
is held at the molecular site by the EP coupling, it is 
trapped to the greatest possible extent at about t/r e ~ 17 
(which almost coincides with the first maximum in (ng) 

at tmax/Te — 20). Self-evidently the maximum is en- 
hanced as e Pi 9 increases. 

When an electron reaches site 9 it can emit a phonon 
to lower its energy. The phonon period is r p h- Hence, 
for the adiabatic regime discussed in Fig. [TJ the phonon 
excited by the first arriving "part of the electron" is still 
present when the last part of the electron arrives (cf. the 
phonon time scale displayed in the graphs at the opposite 
x-axes). Because of this retardation effect the number 
of phonons at the molecular site 9 steadily increases and 
(&9& 9 ) develops, for wo = 0.3, a maximum in time slightly 
before (ng) reaches its maximum. As time proceeds fur- 
ther, the particle starts hopping further away from the 
oscillatory site (recall that the system is no longer in 
an eigenstate after the interaction quench) and (ng) (t) 
decreases until the whole process recurs. Importantly, 
the particle density (ng)(t) at its first mimium around 

*min ~ 2tmax — 40r e is substantially larger than at t = 
above the "critical" EP coupling e Py g/2to ~ 1. This gives 
a first indication that indeed a polaron is formed at site 9, 
in contrast to the weak EP coupling case. 

Figures [H and [3] demonstrate that at larger EP cou- 
plings the particle density at the oscillatory site 9 evolves 
in almost the same manner for a relatively long time 
span. This especially holds for the strong-coupling case 
displayed in Fig. [31 where every incoming electron sticks 
to the molecular site. Thereby the total phonon num- 
ber increases with increasing £ Pj g. Then the interesting 
question is, of course, whether (or to what extent) the 
excited phonons are incorporated in the polaronic quasi- 
particle or rather will be uncorrelated. In our case, where 
the EP coupling acts on a single site only, the electron 
can not carry a phonon cloud away (this-more realistic- 
situation will be investigated in the subsequent two sec- 
tions). Nevertheless we can address this question by an- 
alyzing the phonon distribution function, |c m | 2 (t) with 

Sm=o l c m| 2 W — 1) yielding the weight of the m-phonon 
contribution in the wave function \ip(t))£^ 

Figure 2] gives |c m | 2 (f) for different EP interaction 
strengths, ranging from weak to strong couplings. For 
comparison, the corresponding phonon distribution func- 
tions of the stationary ground states (where e Pi g > Vi) 
is shown. At small £ p ,g, \ip(t)} basically is a zero-phonon 
state at any time. As a matter of course a few phonons 
will be emitted but immediately after will be reabsorbed 
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FIG. 3: (Color online) Time dependence of the particle den- 
sity (upper panel) and phonon number (lower panel) at site 
9 in the strong EP coupling adiabatic regime. Initial state as 
in Fig. [T] £ p ,9 is increased in steps of 0.32. 



as the particle passes the molecular site. Therefore, no 
long-living lattice distortion appears that might trap the 
carrier. The situation dramatically changes as e Pt g ex- 
ceeds the critical coupling strength for polaron forma- 
tion. Now the phonon distribution of the ground state is 
Poisson distributed with maxima at about 4 (e Pi9 = 2), 
9 (e p , 9 = 3.32), and 15 (e p , 9 = 4.92). \tp(t = 20r e )) is a 
multi-phonon state as well. Since ujq is rather small, an 
adiabatic potential (energy) surface emerges that retains 
the incoming electron contributions so that the forma- 
tion of an adiabatic Holstein polaror>2i22 can occur. The 
initial energy of our system (Eq = —1.9696), however, 
does not allow the particle to access the polaronic ground 
state having E (e p , 9 = 2) = -2.521, £ (3.32) = -3.628, 
and £J (4.92) = -5.126. As can be seen from Fig. [H 
the form of the phonon distribution function reflects the 
phonon distribution of excited displaced harmonic oscil- 
lator states, indicating that excited states of the polaron 
were realized instead. These states are known to be sep- 
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FIG. 4: (Color online) Contribution of the m-phonon state 
to \ip(t)} at different times t (given in units r e ). Results for 
various e Pt 9 were compared to the phonon distribution func- 
tion of the system's ground state (e p ,9 > OVi, no interaction 
quench). The phonon frequency is ojq = 0.3. 
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FIG. 5: (Color online) Electron densities on the one- 
dimensional chain at different times t/r e . Initial state as be- 
fore; Wo = 0.3. 



one of the polaron excited states (starting out from the 
well-established polaronic state at e p ,9 = 3.32, see Figs [3] 
and [5]) , the spatiotemporal variation of the various ex- 
cited states is the same, even for very long times and away 
from the central molecular site (cf. the curves marked 
by squares, diamonds and circles for e p ,g = 3.64, 3.96 
and 4.28, respectively). In contrast, if we choose an EP 
coupling which does not match the ground-state energy 
E by lowering the polaron-level ladder, after a while the 
densities evolve quite differently [cf. data for e P:9 = 3.80 
(stars)]. 



arated in energy by lu ~Q Indeed, in going e.g. from 
e P: g = 3.32 to e P: g = 4.92 five additional phonons were 
created (all bound to the polaron), giving rise to a po- 
laron excited state. Note that such kinds of phonon dis- 
tributions were found for the Raman- and infrared-active 
intrinsic localized modes in quasi one-dimensional mixed- 
valence transition-metal complexes^ 

The lower panel of Fig. 0] yields some insight on the 
time scale that the polaron formation process takes place. 
Up to t/T e = 10 the electron radiates successive phonons 
which are still uncorrelated, however. Therefore all 
phonon states below a certain threshold are equally well 
represented in \tp(t)). The zero-phonon state has a larger 
weight, of course, because parts of the electron still re- 
side outside the molecular site. At about t/r e — 15 the 
phonons become correlated, i.e. they are increasingly 
tightly bound to the electron. This process is completed 
at t/r e - 20. 

In Fig. [5] we show snapshots of the electron density 
distribution along the whole chain at various points in 
time. We see that by increasing the EP interaction in 
such a way that the energy of the initial state matches 



2. Anti-adiabatic regime 

We next investigate the limit of large phonon fre- 
quencies. Now T p h is comparable or even smaller than 
r e , which means that a phonon can be excited and re- 
absorbed instantaneously when the electron enters the 
molecular site. During this process the electron will be- 
come (partly) dressed by phonons, provided the EP inter- 
action is sufficiently strong. In this case a non-adiabatic 
Lang-Firsov-type polaron is formed j 33 i 35 This will not 
happen in the weak EP coupling regime illustrated in 
Fig. Because of the extremely large phonon energy, 
only very few phonons can be radiated by the electron 
(see lower panel). The molecular-site particle density 
shown in the upper panel is weakly modulated by the 
phonon emission-absorption processes on the time scale 
Tph and "oscillates" on a time scale related to the system 
size N = 17 (within t/r e = 17 all electrons have visited 
the central site once). No retardation phenomena are ob- 
served (in contrast to Fig.[T]). The dip around t/r e = 10 is 
because electrons continously leaving the dot and those 
arriving at this time mostly come from the boundary 
where the electron density at t = was very small. The 
situation becomes more complex as e Pi g increases and 
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the electron at the molecular site will be partly dressed 
leading to stronger fluctuations of the phonon number. 
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FIG. 6: (Color online) Time dependence of the particle den- 
sity (upper panel) and phonon number (lower panel) at molec- 
ular crystal site 9 in the weak EP coupling anti-adiabatic 
regime. The phonon frequency cjo = 8. The EP coupling is 
switched on at t = 0; different curves belong to e Pl g values 
increased by 0.1. 

At "intermediate" couplings e p ,g = 10 (note that 
e Pi g = lO/wo is of the order of one), the ground-state 
energy Eq(1Q) — —10.1215 and the system comes into 
"resonance" with the initial state (having Eq ~ —2) by 
exciting just one phonon (see Fig. [7]). Since g\ > 1 a po- 
laron will evolve, i.e. the phonon is bound to the electron. 
The same happens at £ Pi9 = 18 (.Bo (18) = -18.0622), 
but now two phonons will be excited and incorporated. 
There are two points worth mentioning. First, as the 
system oscillates on its t/r e — 17-period it can adjust far 
better in order to form a polaron at the sequent t maK - 
points. Second, putting e Pi g only somewhat out of tune, 
both (ng) and (6gfe g ) are substantially reduced for all t, 
i.e. polaron formation is suppressed. 

In the strong EP coupling regime displayed in Fig. [51 
the anti-adiabatic Lang-Firsov polaron has been fully de- 
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FIG. 7: (Color online) Time dependence of the particle 
density (upper panel) and phonon number (lower panel) at 
the molecular crystal site 9 in the strong EP coupling anti- 
adiabatic regime. 



veloped. In this case the arriving electronic contributions 
stay at the molecular site for such a long time that (ng) 

reaches 0.8. Note the increase of £L as compared to 
Figs. [6] and [7] (in order to demonstrate that t£nax is in- 
deed determined by the system size we included results 
for a system with N = 21 sites). The lower panel makes 
clear that now on average many more phonons (oc g|) 
were incorporated. 

The phonon distribution function shown in Fig. [9] cor- 
roborates this scenario. As for the adiabatic case (cf. 
Fig- HI: we observe a transition from an uncorrelated few- 
phonon state to a correlated multi-phonon polaron state. 
The phonon distribution function shows that IV'(tmax)) 
corresponds to an excited polaron with two (four) bound 
phonons at e Pi g = 18 (e Pi g = 33.7). 
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FIG. 8: (Color online) Time dependence of the particle den- 
sity (upper panel) and phonon number (lower panel) at the 
molecular crystal site 9 in the extremely strong EP coupling 
anti-adiabatic regime. The EP interaction is increased in 
steps of cjo- For TV = 21 the particle and phonon numbers 
are displayed at site 11. 



B. Wave-packet injection 

A bare electron injected into a quantum wire coupled 
to the lattice vibrations at every lattice site is another 
instructive example for the quantum dynamics of polaron 
formation.— To this end, at t = 0, we place a Gaussian 
wave packet of width (Tg = 3/ In 10 and momentum K 
centered at site Iq = 4, 



h/,(0)) = A£e [ 



(l-la) 2 



l e iif(*-M c t| 0)) 



(13) 



and let it evolve (A is a normalization constant to ensure 
(^(0)|^(0)) = 1). Again A 4 = Vi. 

Figure [TU] shows snapshots of the local particle densi- 
ties (n,) and phonon numbers (frj&j) for intermediate EP 




FIG. 9: (Color online) Contribution of m-phonon states to 
\ip(t)} at different EP couplings and points of time. Results 
were compared to the stationary phonon distribution function 
of the system's true ground state. The phonon frequency is 
cuo = 8. 




FIG. 10: (Color online) Spatiotemporal evolution of a free 
electron Gaussian wave packet injected with wave vector 
K = n/2 at t = into an 18-site molecular crystal chain with 
PBCs. Model parameters are e p = 1 and uiq — 0.5. Displayed 
is the time evolution (in units of r e ) of the local particle den- 
sities (rii) (black solid lines, open circles), phonon numbers 
{b\bi) (red dot-dashed lines, stars) and local EP correlations 
Xi,i (blue dashed lines). 



couplings (s p = 1 at all sites) and adiabatic phonon fre- 
quencies (u>o = 0.5). In addition we included results for 
the on-site particle-phonon correlations 



Xi,i = ("<(&< + &<■)>• 



(14) 



The particle injected at site 4 is launched to the right 
(K = 7r/2). Shortly after, the electron is not yet dressed 
and moves nearly as fast as a free particle (see black 
curves in the panels for t/r e — 0.02, 1, 2, and 3). 16 At the 
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FIG. 11: (Color online) Spatiotemporal evolution ol a free 
electron Gaussian wave packet injected with wave vector K = 
tt/9 (upper panels) and K — n/2 (lower panels) at t = 0. 
Model parameters are e p — 3 and ljo — 2. Notations as in 
Fig. [10] 



same time the electron emits (creates) phonons along its 
path, in order to reduce its energy to near the bottom of 
the band. In view of the high initial energy Eo(t = 0) = 
and an intermediate EP coupling strength most of the 
phonons radiated are uncorrelated and therefore continue 
to stay near the particle's starting point. Nevertheless 
the particle drags some phonons with it and finally a 
(coherent) polaron wave packet is formed characterized 
by enhanced local particle-phonon correlations (see sites 
5-7 in the panels for t/r e = 6-8; note that the polaronic 
quasiparticle moves with a reduced velocit y 16 ' 25 ) . Owing 
to the moderate EP coupling these signatures are rather 
weak, however, and are further smeared out when the 
polaronic wave packet dissolves in time. 

As Fig. [11] shows, polaron formation becomes more 
pronounced at larger EP couplings (e p = 3; note the 



different scale of the ordinate compared to Fig. fTO)) . even 
if the phonon frequency is enhanced as well (ujq = 2, 
non-adiabatic regime). The phonon distribution and 
enhanced on-site EP correlations indicate that more 
phonons are in the phonon cloud that travels with the 
particle. Thereby the polaron inertial mass is increased. 
Again some unbound phonons stay at the point where the 
particle takes off. While in the upper graph the particle 
is injected with an energy of about the phonon energy 
above the bottom of the band, its energy is much lower 
in the lower graph where K — tt/9. This difference is 
mainly reflected in the number of unbound phonons; in 
both cases the polaronic quasiparticle emerges at about 
t/r e = 8 . . . 10 and shows the same characteristics after- 
wards. 




LU 



e p =3, a) =2 



A £ p =1,to =0.5 


I 




K=7l/2 





.0 10 20 30 40 ,^50 . 




FIG. 12: (Color online) Upper panel: Time dependence of 
the total number of phonons excited in the molecular crystal 
chain as the free electron wave packet evolves into a polaron. 
Lower panel: Time evolution of the kinetic energy part £"kin. 



Since the system is not in an eigenstate we expect to 
find (decaying) oscillations on the time scale of r p h in the 
process of polaron formation; at least if r e and do not 
differ too much and the EP coupling is not too small. 
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This is illustrated by Fig. [121 showing the variation in 
time of the total number of phonons in the system and 
of the kinetic energy part , 



E 



kin 



H.c. 



(15) 



For ujq = 2 (main panels) these oscillations can be 
clearly detected in both quantities. The kinetic energy 
Ekin(£ P — 3,ujq = 2) = —1.316 for the ground state 
of a polaronic system having the same parameters. We 
find that this value can be much better (periodically) 
approached, injecting a particle with lower energy, i.e. 
K = 7r/9. The minima in -Ekin are reached when the 
particle has absorbed some phonons. Afterward the par- 
ticle radiates the phonons again and its kinetic energy in- 
creases. The oscillations are weaker at K = tt/2. While 
the wave vector of the injected wave packet (fT3]) is a 
continuous variable, finite chains with PBC have only a 
finite set of "allowed" K vectors. Because tt/2 is not an 
allowed wave vector of the periodic 18-site system, we in- 
cluded data for K = 4ir/9 (located next to K = tt/2) as 
well, but — as expected — the results do not change qual- 
itatively. The insets give the total phonon number and 
kinetic energy for the parameters of Fig. 1 1 01 i.e. for 
the adiabatic case. Here we can clearly distinguish two 
regimes: until t/r e ~ 6 many unbound phonons were 
created to lower the particle's total energy; then polaron 
formation sets in and the particle attains a kinetic en- 
ergy close to the polaron's ground-state kinetic energy 
^ dn (e p = l,w = 0.5) = -1.823. 



C. Polaron tunneling 

Finally we investigate the tunneling of a polaronic 
quasiparticle through a potential barrier A 12 = A w = 2 
(quantum wall or dot), with additional EP interaction 
e p ,12 = £p,w For i 7^ 12 we fix A$ = 0. The other model 
parameters are choosen to be e p = 1 and ujq = 2. In the 
numerics we account for all states with up to M < 11 
phonons and have checked that in the ground state the 
weight of basis states containing exactly M = 11 phonons 
is, for the largest e P:W , less than 10 -5 . 

The wave packet injected has energy E (0) — and 
moves to the right with K — tt/2 (see Fig.Q2]). We apply 
OBC, so the particle cannot avoid the barrier coming 
from behind. The upper graph describes the situation 
with a barrier at site 12 only. After the polaron is formed 
at t/r e ~ 6 it hits the quantum wall at t/r e ~ 8 — 10 and 
there is mostly reflected. Besides this backscattered par- 
ticle current a minor part of the particle tunnels through 
the barrier, thereby partly stripping and recollecting its 
accompanied phonons (cf . Fig. [T5] below) . Envisaging a 
vibrating molecular quantum dot located at site 12 (mid- 
dle panel), an additive EP interaction e PtW leads to a 
local polaronic level shift that softens the barrier. As 
a result the particle is transmitted to a much greater 
extent than in the former case (compare the results for 
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FIG. 13: (Color online) Quantum dynamics of polaron for- 
mation and polaron tunnelling through a potential barrier 
A12 = A w — 2. Model parameters are e p = 1, uo = 2. The 
barrier or quantum dot is located at site 12 and has a total 
EP coupling (e p + e p , w ). OBCs were used at sites 1 and 18. 
Notation is the same as in Fig. 1111 For further explanation 
see text. 



11 




FIG. 14: (Color online) Time dependence of the particle den- 
sity (upper panel) and phonon number (lower panel) at the 
site of the wall. 




FIG. 15: (Color online) Cumulated particle density to the 
right of the potential wall at different wave vectors K (up- 
per panel) and various EP couplings e p (lower panel), for the 
setup of Fig. [131 



t/r e = 14— 16). If the quantum dot possess a very strong 
EP interaction the polaron digs at the dot site and stays 
there for a long time (see lower panel). Then of course 
both the reflected and transmitted particle current is low. 

Figure [T4l gives the temporal variation of particle den- 
sity and phonon number at the quantum wall/dot site. 
Obviously the phonons somewhat lag behind the elec- 
tron (retardation effect). During the tunneling process 
the phonon number strongly fluctuates. The e PtW = 4 
curve clearly signals the self-trapping of the electron at 
the dot site. The second bump-series is due to electronic 
contributions retaining after being reflected at the sys- 
tem's boundary (site 18). 

The total transmitted electron density is displayed in 
Fig. [TS] for E p>w — A w and various momenta of the in- 
jected wave packet (upper panel) as well as for different 



e p>w at K — tt/2 (lower panel). Of course a higher initial 
energy enhances the transmission through the tunnel bar- 
rier (compare the K = n/2 and K = tt/9 curves). The 
lower panel first of all shows the time delay of the po- 
laron in reaching the barrier compared to a free particle 
(dotted line). More notably, we observe that for e p>w = 2 
the transmission is as high as for free particles, despite 
the fact that the particle is dressed by phonons in a sig- 
nificant way. Note that the dimensionless EP coupling 
parameters at the dot site are 3/2i and g \ 2 = 1.5. This 
points toward the importance of vibration-mediated tun- 
neling processes (doorway vibrons) 



IV. SUMMARY 

In this work, we have presented an efficient numerical 
method to calculate the time evolution of the many-body 
wave function of an interacting electron-phonon system. 
The approach is based on Chebyshev moment expan- 
sion, applied to the time evolution operator. We focused 
on the process of small polaron formation in finite low- 
dimensional quantum structures described by a general- 
ized Holstein Hamiltonian. Both electron and phonon 
quantum dynamics were treated exactly. 

We first started from a non-interacting ground state 
and analyzed the real-time dynamics of the particle den- 
sity and phonon number after a sudden switching-on 
of the electron-phonon coupling at a single oscillatory 
(molecular quantum dot) site. As a consequence of this 
interaction quench the originally free particle can be 
trapped at the "impurity" site after a while. The self- 
trapping process differs in nature for the adiabatic and 
anti-adiabatic regimes of small and large phonon frequen- 
cies, respectively. In the former case, where the phonons 
are slow and retardation effects play an important role, 
a static lattice distortion evolves that causes an effective 
attractive potential for the electron. As a result a Hol- 
stein polaron is formed. In the latter case phonons can 
follow the electron motion almost instantaneously. Hence 
we observe very fast phonon emission and re- absorption 
processes, which — at large EP interaction strengths — 
give rise to a dynamical dressing of the charge carrier 
that enhances the particle's mass and finally leads to its 
immobilization. In both cases the phonon distribution 
function signals the existence of excited bound polaron- 
phonon states. Since our initial state is not an eigenstate 
of the interacting system, we observe the phenomenon of 
recurrence at later times. 

Next we launched a free-electron Gaussian wave packet 
in a one-dimensional system, subjected to EP coupling 
at every site. The injected bare particle is found to ra- 
diate phonons to lower its energy to near the bottom of 
the band. Thereforre part of the phonons stay near the 
electron's starting point and — if the EP coupling is suffi- 
ciently strong — another part of the phonons will be em- 
bedded in a phonon cloud attached to the (moving) par- 
ticle. The latter polaron quasiparticle formation process 
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takes a period of time that depends on the characteris- 
tic electron and phonon times scales, the EP interaction 
strength, and the initial conditions in a very sensitive 
way. We agree with the findings of previous worki£ that 
the question of how long it takes a polaron to form, has 
no simple answer, because there are multiple time scales 
in the dynamics. 

In the last part we investigated the transmission of a 
polaron through a quantum wall or vibrating quantum 
dot. Depending on the barrier height to electron-phonon 
interaction strength ratio, and the characteristic electron 
and phonon times scales, we found opposed behaviors: 
strong reflection; phonon-mediated tunneling; and intrin- 
sic localization of the polaron. Most notably we showed 
that if the polaronic level lowering just compensates the 
repulsive dot potential and the electron and phonon time 
scales are comparable, a rather heavy small polaron, re- 
gardless of its phonon cloud, tunnels like a free electron, 
On the other hand, if there is a mismatch between both 
quantities , we observe strong phonon fluctuations at the 
dot site and transport through the quantum dot becomes 
significantly suppressed. This might motivate further 
investigations of deformable quantum dot systems, e.g. 



with respect to applications as a current switch. 

In conclusion, we have demonstrated that polaron for- 
mation is a subtle non-linear dynamical process which 
is affected by multiple time/energy scales. The proposed 
long-time Chebyshev expansion method — in combination 
with exact diagonalization techniques — is capable of ad- 
dressing such complex problems, which raises the expec- 
tation that our approach can also be used to study the 
time-evolution of quasiparticles in more general situa- 
tions. 
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